#command line arguments
args <- commandArgs(TRUE)

#reading the beagle utility output file of the window based
#data<-read.table("mb1_beagle_util",sep="\t")  #change paramteres
data<-read.table(as.character(args[1]),sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
distance<-read.table("../1000genome/genome1000_4mb1mbedge/1/markers.txt",sep=" ")
colnames(distance)<-c("marker_identifier","dist","allele1","allele2")
merge1<-merge(data,distance,by="marker_identifier")
#reading the beagle utility output file of the whole chr imputation
data<-c()
#data<-read.table("mb4_beagle_util",sep="\t")
data<-read.table(as.character(args[6]),sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
merge3<-merge(data,distance,by="marker_identifier")

#file names
#window<-"1MB"  #change paramteres
window<-as.character(args[2])
allelicr2_name<-paste("whole_default_vs_",window,"_allelic_r_squared.pdf",sep="")
dosager2_name<-paste("whole_default_vs_",window,"_dosage_r_squared.pdf",sep="")
HWE_dosage_r_squared_name<-paste("whole_default_vs_",window,"_HWE_dosage_r_squared.pdf",sep="")

#parameters for this run
window_start<-c(16050035,18050036,20050036,22050036,24050036,26050036,28050036,30050036,32050036,34050036,36050036,38050036,40050036,42050036,44050036,46050036,48050036,50050036)
window_stop<-c(18050036,20050036,22050036,24050036,26050036,28050036,30050036,32050036,34050036,36050036,38050036,40050036,42050036,44050036,46050036,48050036,50050036,51244207)
#actual_window_start<-c(14560202,23560203,33560203,43560203)  #change paramteres
#actual_window_stop<-c(25560203,35560203,45560203,49565873)  #change paramteres
para_start <- as.character(args[4])
para_stop <- as.character(args[5])
actual_window_start<-as.numeric(unlist(strsplit(para_start,',')))
actual_window_stop<-as.numeric(unlist(strsplit(para_stop,',')))





num_windows<-length(window_start)
#edge<-1000000      #change paramteres
edge<-as.numeric(as.character(args[3]))
#removing the edges start & stop of thw window based
merge2<-merge1[which(merge1$dist >=window_start[1]+edge & merge1$dist <= window_stop[num_windows]-edge),]
###actual_window_start[1]<-actual_window_start[1]+edge
###actual_window_stop[num_windows]<-actual_window_stop[num_windows]-edge

#removing the edges start & stop of the whole chr imputation
merge4<-merge3[which(merge3$dist >=window_start[1]+edge & merge3$dist <= window_stop[num_windows]-edge),]

#creating the indicator for different windows of window based
merge2$ind<-0
i<-1
while(i<num_windows+1)
{
	merge2$ind[which(merge2$dist >=window_start[i] & merge2$dist <=window_stop[i] )]<-i	
	i<-i+1
}

#creating the indicator for different windows of whole chr
merge4$ind<-0
i<-1
while(i<num_windows+1)
{
	merge4$ind[which(merge4$dist >=window_start[i] & merge4$dist <=window_stop[i] )]<-i	
	i<-i+1
}


#calculating distance from edges for the window based
nc<-nrow(merge2)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge2$dist[i]-actual_window_start[merge2$ind[i]])
	y<-abs(merge2$dist[i]-actual_window_stop[merge2$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge2$distance<-distance
distance<-c()


#calculating distance from edges for the whole chr based
nc<-nrow(merge4)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge4$dist[i]-actual_window_start[merge4$ind[i]])
	y<-abs(merge4$dist[i]-actual_window_stop[merge4$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge4$distance<-distance
distance<-c()

#plotting allelic r2
pdf(allelicr2_name,width =50,height =50)
plot(merge2$distance, merge2$allelic_r_squared, xlab="Distance from edge", ylab="allelic_r_square", pch='.',col=c("blue"))
points(merge4$distance, merge4$allelic_r_squared, xlab="Distance from edge", ylab="allelic_r_square", pch='.',col=c("red"))
lines(lowess(merge2$distance,merge2$allelic_r_squared), col="black",) 
lines(lowess(merge4$distance,merge4$allelic_r_squared), col="green")
legend("top", horiz=T,col=c("blue","red","black","green"), lty=c(0,0,1,1), c("Window based approach points","4mb window 2mb edge points","Window based approach lowess curve","4mb window 2mb edge lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()

#plotting dosage r2
pdf(dosager2_name,width =50,height =50)
plot(merge2$distance, merge2$dosage_r_squared, xlab="Distance from edge", ylab="dosage_r_square", pch='.',col=c("blue"))
points(merge4$distance, merge4$dosage_r_squared, xlab="Distance from edge", ylab="dosage_r_square", pch='.',col=c("red"))
lines(lowess(merge2$distance,merge2$dosage_r_squared), col="black",) 
lines(lowess(merge4$distance,merge4$dosage_r_squared), col="green")
legend("top", horiz=T,col=c("blue","red","black","green"), lty=c(0,0,1,1), c("Window based approach points","4mb window 2mb edge points","Window based approach lowess curve","4mb window 2mb edge lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()

#plotting HWE dosage r2
pdf(HWE_dosage_r_squared_name,width =50,height =50)
plot(merge2$distance, merge2$HWE_dosage_r_squared, xlab="Distance from edge", ylab="HWE_r_square", pch='.',col=c("blue"))
points(merge4$distance, merge4$HWE_dosage_r_squared, xlab="Distance from edge", ylab="HWE_r_square", pch='.',col=c("red"))
lines(lowess(merge2$distance,merge2$HWE_dosage_r_squared), col="black",) 
lines(lowess(merge4$distance,merge4$HWE_dosage_r_squared), col="green")
legend("top", horiz=T,col=c("blue","red","black","green"), lty=c(0,0,1,1), c("Window based approach points","4mb window 2mb edge points","Window based approach lowess curve","4mb window 2mb edge lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()



